home *** CD-ROM | disk | FTP | other *** search
/ Chip 2007 January, February, March & April / Chip-Cover-CD-2007-02.iso / Pakiet bezpieczenstwa / mini Pentoo LiveCD 2006.1 / mpentoo-2006.1.iso / livecd.squashfs / usr / lib / python2.4 / site-packages / Numeric / MLab.py < prev    next >
Encoding:
Python Source  |  2006-01-20  |  12.0 KB  |  422 lines
  1. """Matlab(tm) compatibility functions.
  2.  
  3. This will hopefully become a complete set of the basic functions available in
  4. matlab.  The syntax is kept as close to the matlab syntax as possible.  One 
  5. fundamental change is that the first index in matlab varies the fastest (as in 
  6. FORTRAN).  That means that it will usually perform reductions over columns, 
  7. whereas with this object the most natural reductions are over rows.  It's perfectly
  8. possible to make this work the way it does in matlab if that's desired.
  9. """
  10. #import Matrix   --- cannot use Matrix module here because Linear Algebra imports it
  11. #                    and Matrix depends on Linear Algebra.
  12.  
  13. from Numeric import *
  14.  
  15. # Elementary Matrices
  16.  
  17. # zeros is from matrixmodule in C
  18. # ones is from Numeric.py
  19.  
  20. import RandomArray
  21. def rand(*args):
  22.     """rand(d1,...,dn) returns a matrix of the given dimensions
  23.     which is initialized to random numbers from a uniform distribution
  24.     in the range [0,1).
  25.     """
  26.     return RandomArray.random(args)
  27.  
  28. def randn(*args):
  29.     """u = randn(d0,d1,...,dn) returns zero-mean, unit-variance Gaussian
  30.     random numbers in an array of size (d0,d1,...,dn)."""
  31.     x1 = RandomArray.random(args)
  32.     x2 = RandomArray.random(args)
  33.     return sqrt(-2*log(x1))*cos(2*pi*x2)
  34.  
  35.  
  36. def eye(N, M=None, k=0, typecode=None):
  37.     """eye(N, M=N, k=0, typecode=None) returns a N-by-M matrix where the 
  38.     k-th diagonal is all ones, and everything else is zeros.
  39.     """
  40.     if M is None: M = N
  41.     if type(M) == type('d'): 
  42.         typecode = M
  43.         M = N
  44.     m = equal(subtract.outer(arange(N), arange(M)),-k)
  45.     return asarray(m,typecode=typecode)
  46.  
  47.  
  48. def tri(N, M=None, k=0, typecode=None):
  49.     """tri(N, M=N, k=0, typecode=None) returns a N-by-M matrix where all
  50.     the diagonals starting from lower left corner up to the k-th are all ones.
  51.     """
  52.     if M is None: M = N
  53.     if type(M) == type('d'): 
  54.         typecode = M
  55.         M = N
  56.     m = greater_equal(subtract.outer(arange(N), arange(M)),-k)
  57.     return m.astype(typecode)
  58.     
  59. # Matrix manipulation
  60.  
  61. def diag(v, k=0):
  62.     """diag(v,k=0) returns the k-th diagonal if v is a matrix or
  63.     returns a matrix with v as the k-th diagonal if v is a vector.
  64.     """
  65.     v = asarray(v)
  66.     s = v.shape
  67.     if len(s)==1:
  68.         n = s[0]+abs(k)
  69.         if k > 0:
  70.             v = concatenate((zeros(k, v.typecode()),v))
  71.         elif k < 0:
  72.             v = concatenate((v,zeros(-k, v.typecode())))
  73.         return eye(n, k=k)*v
  74.     elif len(s)==2:
  75.         v = add.reduce(eye(s[0], s[1], k=k)*v)
  76.         if k > 0: return v[k:]
  77.         elif k < 0: return v[:k]
  78.         else: return v
  79.     else:
  80.             raise ValueError, "Input must be 1- or 2-D."
  81.     
  82.  
  83. def fliplr(m):
  84.     """fliplr(m) returns a 2-D matrix m with the rows preserved and
  85.     columns flipped in the left/right direction.  Only works with 2-D
  86.     arrays.
  87.     """
  88.     m = asarray(m)
  89.     if len(m.shape) != 2:
  90.         raise ValueError, "Input must be 2-D."
  91.     return m[:, ::-1]
  92.  
  93. def flipud(m):
  94.     """flipud(m) returns a 2-D matrix with the columns preserved and
  95.     rows flipped in the up/down direction.  Only works with 2-D arrays.
  96.     """
  97.     m = asarray(m)
  98.     if len(m.shape) != 2:
  99.         raise ValueError, "Input must be 2-D."
  100.     return m[::-1]
  101.     
  102. # reshape(x, m, n) is not used, instead use reshape(x, (m, n))
  103.  
  104. def rot90(m, k=1):
  105.     """rot90(m,k=1) returns the matrix found by rotating m by k*90 degrees
  106.     in the counterclockwise direction.
  107.     """
  108.     m = asarray(m)
  109.     if len(m.shape) != 2:
  110.         raise ValueError, "Input must be 2-D."
  111.     k = k % 4
  112.     if k == 0: return m
  113.     elif k == 1: return transpose(fliplr(m))
  114.     elif k == 2: return fliplr(flipud(m))
  115.     elif k == 3: return fliplr(transpose(m))
  116.  
  117. def tril(m, k=0):
  118.     """tril(m,k=0) returns the elements on and below the k-th diagonal of
  119.     m.  k=0 is the main diagonal, k > 0 is above and k < 0 is below the main
  120.     diagonal.
  121.     """
  122.     svsp = m.spacesaver()
  123.     m = asarray(m,savespace=1)
  124.     out = tri(m.shape[0], m.shape[1], k=k, typecode=m.typecode())*m
  125.     out.savespace(svsp)
  126.     return out
  127.  
  128. def triu(m, k=0):
  129.     """triu(m,k=0) returns the elements on and above the k-th diagonal of
  130.     m.  k=0 is the main diagonal, k > 0 is above and k < 0 is below the main
  131.     diagonal.
  132.     """
  133.     svsp = m.spacesaver()
  134.     m = asarray(m,savespace=1)
  135.     out = (1-tri(m.shape[0], m.shape[1], k-1, m.typecode()))*m
  136.     out.savespace(svsp)
  137.     return out
  138.  
  139. # Data analysis
  140.  
  141. # Basic operations
  142. def max(m,axis=0):
  143.     """max(m,axis=0) returns the maximum of m along dimension axis.
  144.     """
  145.     m = asarray(m)
  146.     return maximum.reduce(m,axis)
  147.  
  148. def min(m,axis=0):
  149.     """min(m,axis=0) returns the minimum of m along dimension axis.
  150.     """
  151.     m = asarray(m)
  152.     return minimum.reduce(m,axis)
  153.  
  154. # Actually from Basis, but it fits in so naturally here...
  155.  
  156. def ptp(m,axis=0):
  157.     """ptp(m,axis=0) returns the maximum - minimum along the the given dimension
  158.     """
  159.     m = asarray(m)
  160.     return max(m,axis)-min(m,axis)
  161.  
  162. def mean(m,axis=0):
  163.     """mean(m,axis=0) returns the mean of m along the given dimension.
  164.        If m is of integer type, returns a floating point answer.
  165.     """
  166.     m = asarray(m)
  167.     return add.reduce(m,axis)/float(m.shape[axis])
  168.  
  169. # sort is done in C but is done row-wise rather than column-wise
  170. def msort(m):
  171.     """msort(m) returns a sort along the first dimension of m as in MATLAB.
  172.     """
  173.     m = asarray(m)
  174.     return transpose(sort(transpose(m)))
  175.  
  176. def median(m):
  177.     """median(m) returns a median of m along the first dimension of m.
  178.     """
  179.     sorted = msort(m)
  180.     if sorted.shape[0] % 2 == 1:
  181.         return sorted[int(sorted.shape[0]/2)]
  182.     else:
  183.         sorted = msort(m)
  184.         index=sorted.shape[0]/2
  185.         return (sorted[index-1]+sorted[index])/2.0
  186.  
  187. def std(m,axis=0):
  188.     """std(m,axis=0) returns the standard deviation along the given 
  189.     dimension of m.  The result is unbiased with division by N-1.
  190.     If m is of integer type returns a floating point answer.
  191.     """
  192.     x = asarray(m)
  193.     n = float(x.shape[axis])
  194.     mx = asarray(mean(x,axis))
  195.     if axis < 0:
  196.         axis = len(x.shape) + axis
  197.     mx.shape = mx.shape[:axis] + (1,) + mx.shape[axis:]
  198.     x = x - mx
  199.     return sqrt(add.reduce(x*x,axis)/(n-1.0))
  200.  
  201. def cumsum(m,axis=0):
  202.     """cumsum(m,axis=0) returns the cumulative sum of the elements along the
  203.     given dimension of m.
  204.     """
  205.     m = asarray(m)
  206.     return add.accumulate(m,axis)
  207.  
  208. def prod(m,axis=0):
  209.     """prod(m,axis=0) returns the product of the elements along the given
  210.     dimension of m.
  211.     """
  212.     m = asarray(m)
  213.     return multiply.reduce(m,axis)
  214.  
  215. def cumprod(m,axis=0):
  216.     """cumprod(m) returns the cumulative product of the elments along the
  217.     given dimension of m.
  218.     """
  219.     m = asarray(m)
  220.     return multiply.accumulate(m,axis)
  221.  
  222. def trapz(y, x=None, axis=-1):
  223.     """trapz(y,x=None,axis=-1) integrates y along the given dimension of
  224.     the data array using the trapezoidal rule.
  225.     """
  226.     y = asarray(y)
  227.     if x is None:
  228.         d = 1.0
  229.     else:
  230.         d = diff(x,axis=axis)
  231.     y = asarray(y)
  232.     nd = len(y.shape)
  233.     slice1 = [slice(None)]*nd
  234.     slice2 = [slice(None)]*nd
  235.     slice1[axis] = slice(1,None)
  236.     slice2[axis] = slice(None,-1)
  237.     return add.reduce(d * (y[slice1]+y[slice2])/2.0,axis)
  238.  
  239. def diff(x, n=1,axis=-1):
  240.     """diff(x,n=1,axis=-1) calculates the n'th difference along the axis specified.
  241.        Note that the result is one shorter in the axis'th dimension.
  242.        Returns x if n == 0. Raises ValueError if n < 0.
  243.     """
  244.     x = asarray(x)
  245.     nd = len(x.shape)
  246.     if nd == 0:
  247.         nd = 1    
  248.     if n < 0:
  249.         raise ValueError, 'MLab.diff, order argument negative.'
  250.     elif n == 0: 
  251.         return x
  252.     elif n == 1:
  253.         slice1 = [slice(None)]*nd
  254.         slice2 = [slice(None)]*nd
  255.         slice1[axis] = slice(1,None)
  256.         slice2[axis] = slice(None,-1)
  257.         return x[slice1]-x[slice2]
  258.     else:
  259.         return diff(diff(x,1,axis), n-1)
  260.  
  261. def cov(m,y=None, rowvar=0, bias=0):
  262.     """Estimate the covariance matrix.
  263.  
  264.     If m is a vector, return the variance.  For matrices where each row
  265.     is an observation, and each column a variable, return the covariance
  266.     matrix.  Note that in this case diag(cov(m)) is a vector of
  267.     variances for each column.
  268.  
  269.     cov(m) is the same as cov(m, m)
  270.  
  271.     Normalization is by (N-1) where N is the number of observations
  272.     (unbiased estimate).  If bias is 1 then normalization is by N.
  273.  
  274.     If rowvar is zero, then each row is a variable with
  275.     observations in the columns.
  276.     """
  277.     if y is None:
  278.         y = m
  279.     else:
  280.         y = y
  281.     if rowvar:
  282.         m = transpose(m)
  283.         y = transpose(y)
  284.     if (m.shape[0] == 1):
  285.         m = transpose(m)
  286.     if (y.shape[0] == 1):
  287.         y = transpose(y)
  288.     N = m.shape[0]
  289.     if (y.shape[0] != N):
  290.         raise ValueError, "x and y must have the same number of observations."    
  291.     m = m - mean(m,axis=0)
  292.     y = y - mean(y,axis=0)
  293.     if bias:
  294.         fact = N*1.0
  295.     else:
  296.         fact = N-1.0
  297.     # 
  298.     val = squeeze(dot(transpose(m),conjugate(y)) / fact)
  299.     return val
  300.  
  301. def corrcoef(x, y=None):
  302.     """The correlation coefficients
  303.     """
  304.     c = cov(x, y)
  305.     d = diag(c)
  306.     return c/sqrt(multiply.outer(d,d))
  307.  
  308.  
  309. # Added functions supplied by Travis Oliphant
  310. import LinearAlgebra
  311. def squeeze(a):
  312.     "squeeze(a) returns a with any ones from the shape of a removed"
  313.     a = asarray(a)
  314.     b = asarray(a.shape)
  315.     return reshape (a, tuple (compress (not_equal (b, 1), b)))
  316.  
  317.  
  318. def kaiser(M,beta):
  319.     """kaiser(M, beta) returns a Kaiser window of length M with shape parameter
  320.     beta. It depends on the cephes module for the modified bessel function i0.
  321.     """
  322.     import cephes
  323.     n = arange(0,M)
  324.     alpha = (M-1)/2.0
  325.     return cephes.i0(beta * sqrt(1-((n-alpha)/alpha)**2.0))/cephes.i0(beta)
  326.  
  327. def blackman(M):
  328.     """blackman(M) returns the M-point Blackman window.
  329.     """
  330.     n = arange(0,M)
  331.     return 0.42-0.5*cos(2.0*pi*n/(M-1)) + 0.08*cos(4.0*pi*n/(M-1))
  332.  
  333.  
  334. def bartlett(M):
  335.     """bartlett(M) returns the M-point Bartlett window.
  336.     """
  337.     n = arange(0,M)
  338.     return where(less_equal(n,(M-1)/2.0),2.0*n/(M-1),2.0-2.0*n/(M-1))
  339.  
  340. def hanning(M):
  341.     """hanning(M) returns the M-point Hanning window.
  342.     """
  343.     n = arange(0,M)
  344.     return 0.5-0.5*cos(2.0*pi*n/(M-1))
  345.  
  346. def hamming(M):
  347.     """hamming(M) returns the M-point Hamming window.
  348.     """
  349.     n = arange(0,M)
  350.     return 0.54-0.46*cos(2.0*pi*n/(M-1))
  351.  
  352. def sinc(x):
  353.     """sinc(x) returns sin(pi*x)/(pi*x) at all points of array x.
  354.     """
  355.     y = pi* where(x == 0, 1.0e-20, x) 
  356.     return sin(y)/y 
  357.  
  358. def eig(v):
  359.     """[x,v] = eig(m) returns the eigenvalues of m in x and the corresponding
  360.     eigenvectors in the rows of v.
  361.     """
  362.     return LinearAlgebra.eigenvectors(v)
  363.  
  364. def svd(v):
  365.     """[u,x,v] = svd(m) return the singular value decomposition of m.
  366.     """
  367.     return LinearAlgebra.singular_value_decomposition(v)
  368.  
  369.  
  370. def angle(z):
  371.     """phi = angle(z) return the angle of complex argument z."""
  372.     z = asarray(z)
  373.     if z.typecode() in ['D','F']:
  374.        zimag = z.imag
  375.        zreal = z.real
  376.     else:
  377.        zimag = 0
  378.        zreal = z
  379.     return arctan2(zimag,zreal)
  380.  
  381.  
  382. def roots(p):
  383.     """ return the roots of the polynomial coefficients in p.
  384.  
  385.     The values in the rank-1 array p are coefficients of a polynomial.
  386.     If the length of p is n+1 then the polynomial is
  387.     p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]
  388.     """
  389.     if type(p) in [types.IntType, types.FloatType, types.ComplexType]:
  390.         p = asarray([p])
  391.     else:
  392.         p = asarray(p)
  393.  
  394.     n = len(p)
  395.     if len(p.shape) != 1:
  396.         raise ValueError, "Input must be a rank-1 array."
  397.  
  398.     # Strip zeros at front of array
  399.     ind = 0
  400.     while (p[ind] == 0):
  401.         ind = ind + 1
  402.     p = asarray(p[ind:])
  403.  
  404.     N = len(p)
  405.     root = zeros((N-1,),'D')
  406.     # Strip zeros at end of array which correspond to zero-valued roots
  407.     ind = len(p)
  408.     while (p[ind-1]==0):
  409.         ind = ind - 1
  410.     p = asarray(p[:ind])
  411.  
  412.     N = len(p)
  413.     if N > 1:
  414.         A = diag(ones((N-2,),p.typecode()),-1)
  415.         A[0,:] = -p[1:] / p[0]
  416.         root[:N-1] = eig(A)[0]
  417.     if ((root.typecode() == 'F' and allclose(root.imag, 0, rtol=1e-7)) or
  418.         (root.typecode() == 'D' and allclose(root.imag, 0, rtol=1e-14))):
  419.         root = root.real
  420.     return root
  421.  
  422.